#include <iostream>
#include <cstring>
#include <cmath>
#include <algorithm>
#include <utility>
#include <map>

#include "TH1D.h"
#include "TH2D.h"
#include "TProfile.h"
#include "TFile.h"
#include "TTree.h"
#include "TChain.h"
#include "TGraph.h"
#include "TParameter.h"
#include "DelphesClasses.h"
#include "ExRootTreeReader.h"
#include "mr.h"
#include "TClonesArray.h"
#include "covariance.h"
#include "vector_tools.h"

#include "CmdLine/CmdLine.hh"
using namespace std;

bool debug=false;

const double PI=3.1415926;

class TLorentzSorter{
public:
  bool operator()(const TLorentzVector& a, const TLorentzVector& b) const
  {
    return a.Pt()>b.Pt();
  }
};


double fixphi(double input)
{
  if(input < 0)
    input = -input;

  while(input > 2*PI)
    {
      input -= 2*PI;
    }
  return input;
}

//send deltaphi = PI to infinity
double transformphi(double input)
{
  return log(tan(fixphi(input)/2));
}


//declare a tree class
class MyTree{
public:
   TChain         *t;   //!pointer to the analyzed TTree or TChain

   // Declaration of leaf types
  double weight_scale;

   Float_t         weight;
   Float_t         MET;
   Float_t         MET_phi;
   Float_t         HTjets;
   Float_t         HTleptons;
   Float_t         HTleptons_lep20;
   Float_t         m_eff;
   Float_t         METsignif;
   Int_t           nBtags;
   Int_t           nJets30;
   Int_t           nJets60;
   Int_t           nJets100;
   Int_t           nJets200;
   Float_t         pTJet[100];
   Float_t         etaJet[100];
   Float_t         phiJet[100];
   Float_t         massJet[100];
   UInt_t          btagJet[100];

   Float_t         Tau1Jet[100];
   Float_t         Tau2Jet[100];
   Float_t         Tau3Jet[100];
   Int_t           WTagJet[100];
   Float_t         TrimmedMassJet[100];
   Int_t           TopTagJet[100];

   Float_t         mindPhiMetJet30;
   Float_t         mindPhiMetJet60;
   Float_t         mindPhiMetJet100;
   Float_t         mindPhiMetJet200;
   Float_t         dPhiMETBtag0;
   Float_t         dPhiMETBtag1;
   Float_t         dPhiMETBtag2;
   Float_t         mindPhiMETBtag;
   Int_t           nLeps35;
   Float_t         pTLep[10];
   Float_t         etaLep[10];
   Float_t         phiLep[10];
   Int_t           flavLep[10];
   Float_t         pTTruth[50];
   Float_t         etaTruth[100];
   Float_t         phiTruth[100];
   Float_t         massTruth[100];
   Int_t           pdgIdTruth[100];
   Int_t           M1Truth[100];
   Int_t           M2Truth[100];
   Int_t           D1Truth[100];
   Int_t           D2Truth[100];
   Int_t           statusTruth[100];
  


  MyTree(string input)
  {
    // Set branch addresses and branch pointers
    t = new TChain("OutputTree");

    t->AddFile(input.c_str());
    assert(t->GetFile() != NULL);

    t->SetMakeClass(1);
    
    t->SetBranchAddress("weight", &weight);
    t->SetBranchAddress("MET", &MET);
    t->SetBranchAddress("MET_phi", &MET_phi);
    t->SetBranchAddress("HTjets", &HTjets);
    t->SetBranchAddress("HTleptons", &HTleptons);
    t->SetBranchAddress("HTleptons_lep20", &HTleptons_lep20);
    t->SetBranchAddress("m_eff", &m_eff);
    t->SetBranchAddress("METsignif", &METsignif);
    t->SetBranchAddress("nBtags", &nBtags);
    t->SetBranchAddress("nJets30", &nJets30);
    t->SetBranchAddress("nJets60", &nJets60);
    t->SetBranchAddress("nJets100", &nJets100);
    t->SetBranchAddress("nJets200", &nJets200);
    t->SetBranchAddress("pTJet", pTJet);
    t->SetBranchAddress("etaJet", etaJet);
    t->SetBranchAddress("phiJet", phiJet);
    t->SetBranchAddress("massJet", massJet);
    t->SetBranchAddress("btagJet", btagJet);
    t->SetBranchAddress("Tau1Jet", Tau1Jet);
    t->SetBranchAddress("Tau2Jet", Tau2Jet);
    t->SetBranchAddress("Tau3Jet", Tau3Jet);
    t->SetBranchAddress("WTagJet", WTagJet);
    t->SetBranchAddress("TrimmedMassJet", TrimmedMassJet);
    t->SetBranchAddress("TopTagJet", TopTagJet);

    t->SetBranchAddress("mindPhiMetJet30", &mindPhiMetJet30);
    t->SetBranchAddress("mindPhiMetJet60", &mindPhiMetJet60);
    t->SetBranchAddress("mindPhiMetJet100", &mindPhiMetJet100);
    t->SetBranchAddress("mindPhiMetJet200", &mindPhiMetJet200);
    t->SetBranchAddress("dPhiMETBtag0", &dPhiMETBtag0);
    t->SetBranchAddress("dPhiMETBtag1", &dPhiMETBtag1);
    t->SetBranchAddress("dPhiMETBtag2", &dPhiMETBtag2);
    t->SetBranchAddress("mindPhiMETBtag", &mindPhiMETBtag);
    t->SetBranchAddress("nLeps35", &nLeps35);
    t->SetBranchAddress("pTLep", pTLep);
    t->SetBranchAddress("etaLep", etaLep);
    t->SetBranchAddress("phiLep", phiLep);
    t->SetBranchAddress("flavLep", flavLep);
    t->SetBranchAddress("pTTruth", pTTruth);
    t->SetBranchAddress("etaTruth", etaTruth);
    t->SetBranchAddress("phiTruth", phiTruth);
    t->SetBranchAddress("massTruth", massTruth);
    t->SetBranchAddress("pdgIdTruth", pdgIdTruth);
    t->SetBranchAddress("M1Truth", M1Truth);
    t->SetBranchAddress("M2Truth", M2Truth);
    t->SetBranchAddress("D1Truth", D1Truth);
    t->SetBranchAddress("D2Truth", D2Truth);
    t->SetBranchAddress("statusTruth", statusTruth);


    //now if we have signal, need to change
    weight_scale=-1.0;
    //this is signal file
    if(input.find("decoupledGO")!=string::npos)
      {
	//get the cxn
	TParameter<double>* par=
	  (TParameter<double>*) t->GetFile()->Get("TotalProcessEvents");
	//divide by the total number of generated events
	double nevt=par->GetVal();


	//get a TGraph for the cxn
	double mass[]={50,200,400,600,800,1000,1200,1400,1600,1800,2000,2200,2400,2600,3400,3600,3800,4000,4200,4400,4600,4800,5000,5200,5400,5600,5800,6000,6200,6400,6600,6800,7000,7200,7400,7600,7800,8000,8200,8400,8600,8800,9000,9200,9400,9600,9800,10000, 10200, 10400, 10600, 10800, 11000, 11200, 11400, 11600, 11800, 12000, 12200, 12400, 12600, 12800, 13000, 13200, 13400, 13600, 13800, 14000, 14200, 14400, 14600, 14800, 15000, 15200, 15400, 15600, 15800, 16000, 16200, 16400, 16600, 16800, 17000, 17200, 17400, 17600, 17800, 18000, 18200, 18400, 18600, 18800, 19000, 19200, 19400, 19600, 19800, 20000};


	double cxn[]={4.180000e+05, 3.000000e+03, 1.950000e+02, 3.660000e+01, 1.020000e+01, 3.960000e+00, 1.810000e+00, 9.450000e-01, 5.210000e-01, 2.910000e-01, 1.794000e-01, 1.090000e-01, 6.860000e-02, 4.440000e-02, 9.880000e-03, 7.100000e-03, 5.160000e-03, 3.800000e-03, 2.840000e-03, 2.140000e-03, 1.622000e-03, 1.244000e-03, 9.600000e-04, 7.460000e-04, 5.840000e-04, 4.600000e-04, 3.660000e-04, 2.900000e-04, 2.320000e-04, 1.872000e-04, 1.512000e-04, 1.226000e-04, 9.980000e-05, 8.140000e-05, 6.680000e-05, 5.480000e-05, 4.520000e-05, 3.740000e-05, 3.100000e-05, 2.580000e-05, 2.140000e-05, 1.788000e-05, 1.496000e-05, 1.254000e-05, 1.052000e-05, 8.840000e-06, 7.460000e-06, 6.280000e-06, 5.340000e-06, 4.520000e-06, 3.860000e-06, 3.280000e-06, 2.800000e-06, 2.380000e-06, 2.040000e-06, 1.742000e-06, 1.492000e-06, 1.278000e-06, 1.096000e-06, 9.400000e-07, 8.080000e-07, 6.940000e-07, 5.980000e-07, 5.140000e-07, 4.420000e-07, 3.820000e-07, 3.300000e-07, 2.840000e-07, 2.440000e-07, 2.120000e-07, 1.828000e-07, 1.580000e-07, 1.366000e-07, 1.182000e-07, 1.022000e-07, 8.840000e-08, 7.660000e-08, 6.620000e-08, 5.740000e-08, 4.980000e-08, 4.300000e-08, 3.740000e-08, 3.240000e-08, 2.800000e-08, 2.440000e-08, 2.100000e-08, 1.828000e-08, 1.584000e-08, 1.374000e-08, 1.190000e-08, 1.032000e-08, 8.940000e-09, 7.760000e-09, 6.720000e-09, 5.820000e-09, 5.040000e-09, 4.380000e-09, 3.800000e-09};




	TGraph g(sizeof(mass)/sizeof(double), mass, cxn);
	
	//now get the mass
	int pos=input.find("Up.T1");
	
	string sub_mass=input.substr(pos+5, 5);
	cout<<"found mass string for stop production: "<<sub_mass<<endl;

	int real_mass = atoi(sub_mass.c_str());

	double real_cxn = g.Eval(real_mass);
	cout<<"cxn: "<<real_cxn<<endl;

	weight_scale= real_cxn * 1.0e3 / nevt;

	
	cout<<"squark mass: "<<real_mass<<endl;
	cout<<"weight: "<<weight_scale<<endl;
      }

  }


  int cout_truth() const
  {
    for(int i=0; i<50; i++)
      cout<<i<<": \t"<<pdgIdTruth[i]
	  <<"\t"<<M1Truth[i]-2
	  <<"\t"<<statusTruth[i]<<endl;
  }

  //function to find parent of muon
  int mu_parent(int i) const
  {
    if(i < 0 || i >= 50)
      return -1;

    int parent=M1Truth[i]-2;
    
    while(parent >= 0 && parent <50 && 
	  abs(pdgIdTruth[parent])==13 )
      {
	int temp_parent=parent;
	parent=M1Truth[parent]-2; 
	
	//infinite loop
	if(temp_parent == parent)
	  break;
      }

    if(parent >= 0 && parent <50)
      return abs(pdgIdTruth[parent]);

    else
      return -1;
  }

  //supply a vector denoting the charge
  vector<TLorentzVector> get_muons(double pt_cut=100) const
  {
    vector<TLorentzVector> mu_list;
    
    for(int ic=0; ic<50*2; ++ic)
      {
	int i=ic;
	if(ic>=50) 
	  i-=50;
	
	if(abs(pdgIdTruth[i])!=13) continue;
	if(pTTruth[i] < pt_cut) continue;
	if(fabs(etaTruth[i])>2.5) continue;


	//do status 2 muons in second loop
	if(ic < 50 && statusTruth[i]==2) 
	  continue;

	TLorentzVector muj;

	muj.SetPtEtaPhiM(pTTruth[i], etaTruth[i], 
			 phiTruth[i], 0);
	
	
	//check for duplicates
	bool b_dup=false;
	for(int k=0; k<mu_list.size(); ++k)
	  {
	    if(mu_list[k].DeltaR(muj) < 1e-3)
	      {
		b_dup=true;
		break;
	      }
	  }
	if(b_dup)
	  continue;

	mu_list.push_back(muj);
      }

    //sort the muons from decreasing pt
    sort(mu_list.begin(), mu_list.end(), TLorentzSorter());

    return mu_list;
  }


  //return all the top quarks
  vector<TLorentzVector> get_particle(int pdgid=6) const
  {
    vector<TLorentzVector> top_list;
    
    for(int i=0; i<50; ++i)
      {
	
	if(abs(pdgIdTruth[i])!=pdgid) continue;

	TLorentzVector topj;

	topj.SetPtEtaPhiM(pTTruth[i], etaTruth[i], 
			 phiTruth[i], 0);
	
	
	top_list.push_back(topj);
      }

    //sort the muons from decreasing pt
    sort(top_list.begin(), top_list.end(), TLorentzSorter());

    return top_list;
  }


};


bool match_jet(const TLorentzVector& jet, 
	       const vector<TLorentzVector>& particles,
	       double delta_R=0.5)
{

  bool good_match=false;
  for(int i=0; i<particles.size(); i++)
    {
      if(particles[i].DeltaR(jet) < delta_R)
	{
	  good_match=true;
	  break;
	}
    }
  return good_match;
}


int main(int argc, char** argv)
{
  //initialize input
  CmdLine cmdline(argc, argv);

  //get input and output
  string input, output;

  double MR_cut, met_cut;

  bool b_W_only=false, b_Z_only=false;
  
  try
    {
      input = cmdline.value<string> ("-i");
      output = cmdline.value<string> ("-o");
      MR_cut = cmdline.value<double> ("-MR", 6000);
      met_cut = cmdline.value<double> ("-met", 4000);

      b_W_only = cmdline.present ("-W");
      b_Z_only = cmdline.present ("-Z");

    }
  catch(...)
    {
      cout<<"ERROR: please specify input and output via -i input -o output"<<endl;
      return 0;
    }

  bool b_compressed = met_cut <10;

  //initialize plots
  map<string, TH1D*> hist;
  map<string, TH2D*> hist2D;

  //book histograms

  hist["cutflow"]=new TH1D("cutflow", "Cut Flow", 10, 0, 10);

  //cutflow 2 = remove jet mass cut
  hist["cutflow2"]=new TH1D("cutflow2", "Cut Flow 2", 10, 0, 10);
  //cutflow 3 = remove jet mass cut + raise pt cut to 1000
  hist["cutflow3"]=new TH1D("cutflow3", "Cut Flow 3", 10, 0, 10);
  //cutflow 4 = remove jet mass cut + raise pt cut to 1000
  // + change to met cut of 3000
  hist["cutflow4"]=new TH1D("cutflow4", "Cut Flow 4", 10, 0, 10);
  //cutflow 5 = same as cutflow without jet mass, but with 2 muons
  hist["cutflow5"]=new TH1D("cutflow5", "Cut Flow 5", 10, 0, 10);

  //cutflow 5 = same as cutflow but with muon pt cut > 500 and jet pt 1000
  hist["cutflow6"]=new TH1D("cutflow6", "Cut Flow 6", 10, 0, 10);

  hist["cutflow_old"]=new TH1D("cutflow_old", "Cut Flow", 10, 0, 10);
  hist["fisher"]=new TH1D("fisher", "Fisher Discriminant", 200, -5, 5);

  hist["mu1_parent"]=new TH1D("mu1_parent", "first #mu parent", 10, 0, 10);

  hist["mu2_parent"]=new TH1D("mu2_parent", "second #mu parent", 10, 0, 10);


  hist["nj100"]=new TH1D("nj100", "Number of Jets with p_{T} > 100 GeV", 20, 0, 20);

  hist["nmuj"]=new TH1D("nmuj", "Number of Muons within jets", 20, 0, 20);
  hist["nmu"]=new TH1D("nmu", "Number of Muons", 20, 0, 20);

  hist["nmu_gen"]=new TH1D("nmu_gen", "Number of Muons", 20, 0, 20);

  hist["nW_gen"]=new TH1D("nW_gen", "Number of W", 20, 0, 20);
  hist["nZ_gen"]=new TH1D("nZ_gen", "Number of Z", 20, 0, 20);


  hist["nmu_gen"]=new TH1D("nmu_gen", "Number of mu", 20, 0, 20);
  hist["ne_gen"]=new TH1D("ne_gen", "Number of e", 20, 0, 20);
  hist["ntau_gen"]=new TH1D("ntau_gen", "Number of tau", 20, 0, 20);
  hist["nneu_gen"]=new TH1D("nneu_gen", "Number of neu", 20, 0, 20);

  hist["nmu_after"]=new TH1D("nmu_after", "Number of mu", 20, 0, 20);
  hist["ne_after"]=new TH1D("ne_after", "Number of e", 20, 0, 20);
  hist["ntau_after"]=new TH1D("ntau_after", "Number of tau", 20, 0, 20);
  hist["nneu_after"]=new TH1D("nneu_after", "Number of neu", 20, 0, 20);



  hist["mu_parent"]=new TH1D("nmu_parent", "Number of Muons (parents)", 20, 0, 20);

  hist["nj100_after"]=new TH1D("nj100_after", "Number of Jets with p_{T} > 100 GeV", 20, 0, 20);

  hist["nbj100"]=new TH1D("nbj100", "Number of B Jets with p_{T} > 100 GeV", 10, 0, 10);
  hist["nbj100_after"]=new TH1D("nbj100_after", "Number of B Jets with p_{T} > 100 GeV", 20, 0, 20);


  hist["nj150"]=new TH1D("nj150", "Number of Jets with p_{T} > 150 GeV", 20, 0, 20);

  hist["HT"]=new TH1D("HT", "H_{T}", 20000/200, 0, 20000);
  hist["HT_after"]=new TH1D("HT_after", "H_{T} (after cut)", 20000/200, 0, 20000);

  hist["MR"]=new TH1D("MR", "M_{R}", 20000/200, 0, 20000);
  hist["MR_after"]=new TH1D("MR_after", "M_{R} (after cut)", 20000/200, 0, 20000);

  hist["MJ"]=new TH1D("MJ", "M_{J}", 20000/200, 0, 20000);

  hist["met"]=new TH1D("met", "#slash{E}_{T}", 10000/200, 0, 10000);
  hist["metsig"]=new TH1D("metsig", "#slash{E}_{T}/#sqrt{H_{T}}", 200, 0, 150);

  hist["met_after"]=new TH1D("met_after", "#slash{E}_{T}", 10000/200, 0, 10000);
  hist["metsig_after"]=new TH1D("metsig_after", "#slash{E}_{T}/#sqrt{H_{T}} (after cut)", 200, 0, 150);


  hist["jet1pt"]=new TH1D("jet1pt", "1^{st} jet p_{T}", 10000/200, 0, 10000);
  hist["jet1pt_after"]=new TH1D("jet1pt_after", "1^{st} jet p_{T} (after cut)", 10000/200, 0, 10000);
  hist["jet2pt"]=new TH1D("jet2pt", "2^{nd} jet p_{T}", 10000/200, 0, 10000);
  hist["jet2pt_after"]=new TH1D("jet2pt_after", "2^{nd} jet p_{T} (after cut)", 10000/200, 0, 10000);
  hist["jet3pt"]=new TH1D("jet3pt", "3^{rd} jet p_{T}", 10000/200, 0, 10000);
  hist["jet4pt"]=new TH1D("jet4pt", "4^{th} jet p_{T}", 10000/200, 0, 10000);


  hist["jet1mass"]=new TH1D("jet1mass", "1^{st} jet mass", 1000/20, 0, 1000);
  hist["jet1mass_after"]=new TH1D("jet1mass_after", "1^{st} jet mass (after cut)", 1000/20, 0, 1000);
  hist["jet2mass"]=new TH1D("jet2mass", "2^{nd} jet mass", 1000/20, 0, 1000);
  hist["jet2mass_after"]=new TH1D("jet2mass_after", "2^{nd} jet mass (after cut)", 1000/20, 0, 1000);
  hist["jet3mass"]=new TH1D("jet3mass", "3^{rd} jet mass", 1000/20, 0, 1000);
  hist["jet4mass"]=new TH1D("jet4mass", "4^{th} jet mass", 1000/20, 0, 1000);

  
  hist["dphi"]=new TH1D("dphi", "min #Delta #phi between #slash{E}_{T} and all jet", 100, 0, PI);

  hist["dphi_after"]=new TH1D("dphi_after", "min #Delta #phi between #slash{E}_{T} and all jet (after cut)", 100, 0, PI);


  //a bunch of distributions after some mild cut
  hist["HT_mild"]=new TH1D("HT_mild", "H_{T}", 20000/200, 0, 20000);
  hist["MR_mild"]=new TH1D("MR_mild", "M_{R}", 20000/200, 0, 20000);
  hist["MRHT_mild"]=new TH1D("MRHT_mild", "M_{R}/H_{T}", 100, 0, 10);
  
  hist["nj100_mild"]=new TH1D("nj100_mild", "Number of Jets with p_{T} > 100 GeV", 20, 0, 20);

  hist["nl35_mild"]=new TH1D("nl35_mild", "Number of Leptons with p_{T} > 35 GeV", 20, 0, 20);

  hist["ptmuj_mild"]=new TH1D("ptmuj_mild", "#mu/jet p_{T} ratio (mild cut)", 50, 0, 1.0);

  hist["dphimin_mild"]=new TH1D("dphimin_mild", "min #Delta #phi between #slash{E}_{T} and jets", 50, 0, PI);

  hist["dphimin500_mild"]=new TH1D("dphimin500_mild", "min #Delta #phi between #slash{E}_{T} and 500 GeV jets", 50, 0, PI);

  hist["dphimin200_mild"]=new TH1D("dphimin200_mild", "min #Delta #phi between #slash{E}_{T} and 200 GeV jets", 50, 0, PI);

  hist["dphi1_mild"]=new TH1D("dphi1_mild", "#Delta #phi between #slash{E}_{T} and 1^{st} jet", 50, 0, PI);

  hist["dphi2_mild"]=new TH1D("dphi2_mild", "#Delta #phi between #slash{E}_{T} and 2^{nd} jet", 50, 0, PI);

  hist["dphimin_mod_mild"]=new TH1D("dphimin_mod_mild", "min #Delta #phi between #slash{E}_{T} and jets", 200, -10, 10);

  hist["dphimin500_mod_mild"]=new TH1D("dphimin500_mod_mild", "min #Delta #phi between #slash{E}_{T} and 500 GeV jets", 200, -10, 10);

  hist["dphimin200_mod_mild"]=new TH1D("dphimin200_mod_mild", "min #Delta #phi between #slash{E}_{T} and 200 GeV jets", 200, -10, 10);

  hist["dphi1_mod_mild"]=new TH1D("dphi1_mod_mild", "#Delta #phi between #slash{E}_{T} and 1^{st} jet", 200, -10, 10);

  hist["dphi2_mod_mild"]=new TH1D("dphi2_mod_mild", "#Delta #phi between #slash{E}_{T} and 2^{nd} jet", 200, -10, 10);

  hist["dR_mild"]=new TH1D("dR_mild", "#DeltaR between 1^{st} and 2^{nd} jet", 50, 0, 2*PI);

  hist["dphij_mild"]=new TH1D("dphij_mild", "#Delta#phi between 1^{st} and 2^{nd} jet", 50, 0, PI);
  hist["deta_mild"]=new TH1D("deta_mild", "#Delta#eta between 1^{st} and 2^{nd} jet", 50, 0, 5);
  
  hist["jet1pt_mild"]=new TH1D("jet1pt_mild", "1^{st} jet p_{T} (mild cut)", 10000/200, 0, 10000);

  hist["mjj_mild"]=new TH1D("mjj_mild", "Leading dijet mass", 10000/200, 0, 10000);

  hist["jet2pt_mild"]=new TH1D("jet2pt_mild", "2^{nd} jet p_{T} (mild cut)", 10000/200, 0, 10000);

  hist["jetptratio_mild"]=new TH1D("jetptratio_mild", "jet p_{T} ratio (mild cut)", 50, 0, 1.0);

  hist["jet1mass_mild"]=new TH1D("jet1mass_mild", "1^{st} jet mass (mild cut)", 1000/20, 0, 1000);

  hist["jet2mass_mild"]=new TH1D("jet2mass_mild", "2^{nd} jet mass (mild cut)", 1000/20, 0, 1000);

  hist["jet1masstrim_mild"]=new TH1D("jet1masstrim_mild", "1^{st} jet trimmed mass (mild cut)", 1000/20, 0, 1000);

  hist["jet2masstrim_mild"]=new TH1D("jet2masstrim_mild", "2^{nd} jet trimmed mass (mild cut)", 1000/20, 0, 1000);

  hist["met_mild"]=new TH1D("met_mild", "#slash{E}_{T}", 10000/100, 0, 10000);

  hist["metsig_mild"]=new TH1D("metsig_mild", "#slash{E}_{T}/#sqrt{H_{T}}", 200, 0, 150);

  hist["nbj100_mild"]=new TH1D("nbj100_mild", "Number of B Jets with p_{T} > 100 GeV", 10, 0, 10);

  hist["jet1t32_mild"]=new TH1D("jet1t32_mild", "1^{st} jet #tau_{32}", 25, 0, 1);

  hist["jet2t32_mild"]=new TH1D("jet2t32_mild", "2^{nd} jet #tau_{32}", 25, 0, 1);

  hist["jet1t31_mild"]=new TH1D("jet1t31_mild", "1^{st} jet #tau_{31}", 25, 0, 1);

  hist["jet2t31_mild"]=new TH1D("jet2t31_mild", "2^{nd} jet #tau_{31}", 25, 0, 1);


  //allbut plots
  hist["jet1pt_allbut"]=new TH1D("jet1pt_allbut", "1^{st} jet p_{T} (all but)", 10000/200, 0, 10000);

  hist["jet2pt_allbut"]=new TH1D("jet2pt_allbut", "2^{nd} jet p_{T} (all but)", 10000/200, 0, 10000);

  hist["HT_allbut"]=new TH1D("HT_allbut", "H_{T} (all but)", 20000/200, 0, 20000);

  hist["nmuj_allbut"]=new TH1D("nmuj_allbut", "Number of Muons within jets", 20, 0, 20);

  hist["nl35_allbut"]=new TH1D("nl35_allbut", "Number of isolated leptons", 21, -1, 20);

  hist["nl100_allbut"]=new TH1D("nl100_allbut", "Number of isolated leptons", 20, 0, 20);

  hist["nl200_allbut"]=new TH1D("nl200_allbut", "Number of isolated leptons", 20, 0, 20);

  hist["nmuj500_allbut"]=new TH1D("nmuj500_allbut", "Number of Muons within jets", 20, 0, 20);

  hist["nmuj200_allbut"]=new TH1D("nmuj200_allbut", "Number of Muons within jets", 20, 0, 20);

  hist["jet1masstrim_allbut"]=new TH1D("jet1masstrim_allbut", "1^{st} jet trimmed mass (all but)", 1000/20, 0, 1000);

  hist["jet2masstrim_allbut"]=new TH1D("jet2masstrim_allbut", "2^{nd} jet trimmed mass (all but)", 1000/20, 0, 1000);

  hist["met_allbut"]=new TH1D("met_allbut", "#slash{E}_{T}", 10000/100, 0, 10000);

  hist["met_allbut_dilep_dR"]=new TH1D("met_allbut_dilep_dR", "#slash{E}_{T}", 10000/100, 0, 10000);

  hist["met_allbut_dilep_dR_loose"]=new TH1D("met_allbut_dilep_dR_loose", "#slash{E}_{T}", 10000/100, 0, 10000);

  hist["met_allbut_dilep_mod"]=new TH1D("met_allbut_dilep_mod", "#slash{E}_{T}", 10000/100, 0, 10000);


  hist["met_allbut_dilep"]=new TH1D("met_allbut_dilep", "#slash{E}_{T}", 10000/100, 0, 10000);
  hist["met_allbut_dilep_nolep"]=new TH1D("met_allbut_dilep_nolep", "#slash{E}_{T}", 10000/100, 0, 10000);

  hist["met_allbut_mild"]=new TH1D("met_allbut_mild", "#slash{E}_{T}", 10000/100, 0, 10000);
  hist["met_allbut_milder"]=new TH1D("met_allbut_milder", "#slash{E}_{T}", 10000/100, 0, 10000);

  hist["met_allbut_nolep"]=new TH1D("met_allbut_nolep", "#slash{E}_{T}", 10000/100, 0, 10000);

  hist["metsig_allbut"]=new TH1D("metsig_allbut", "#slash{E}_{T}/#sqrt{H_{T}}", 200, 0, 150);

  hist["dphimin200_allbut"]=new TH1D("dphimin200_allbut", "min #Delta #phi between #slash{E}_{T} and 200 GeV jets", 50, 0, PI);

  hist["dphilep_allbut"]=new TH1D("dphilep_allbut", "min #Delta #phi between #slash{E}_{T} and leptons", 50, 0, PI);

  hist["dphimin500_allbut"]=new TH1D("dphimin500_allbut", "min #Delta #phi between #slash{E}_{T} and 500 GeV jets", 50, 0, PI);


  hist["mjj_after"]=new TH1D("mjj_after", "Leading dijet mass", 10000/200, 0, 10000);

  hist["mll_after"]=new TH1D("mll_after", "dilepton mass", 1000/10, 0, 1000);

  hist["ptl1_after"]=new TH1D("ptl1_after", "Leading Lepton p_{T}", 1000/10, 0, 1000);

  hist["ptl2_after"]=new TH1D("ptl2_after", "Second Lepton p_{T}", 1000/10, 0, 1000);

  hist["dphill_after"]=new TH1D("dphill_after", "lepton d#phi", 100, 0, PI);
  
  hist["dRll_after"]=new TH1D("dRll_after", "lepton dR", 200, 0, 4);

  hist["ST_after"]=new TH1D("ST_after", "sum of all p_{T}", 20000/200, 0, 20000);


  hist["top_pt1"]=new TH1D("top_pt1", "Leading top p_{T}", 10000/100, 0, 10000);

  hist["top_pt2"]=new TH1D("top_pt2", "Second top p_{T}", 10000/100, 0, 10000);

  hist["muon_pt1"]=new TH1D("muon_pt1", "Leading muon p_{T}", 10000/100, 0, 10000);


  hist["top_dR"]=new TH1D("top_dR", "top #DeltaR", 100, 0, 5);
  hist["top_dR_after"]=new TH1D("top_dR_after", "top #DeltaR (after cut)", 100, 0, 5);

  hist["top_pt1_after"]=new TH1D("top_pt1_after", "Leading top p_{T}", 10000/100, 0, 10000);

  hist["top_pt2_after"]=new TH1D("top_pt2_after", "Second top p_{T}", 10000/100, 0, 10000);



  //initialize profiles
  map<string, TProfile*> prof;


  prof["eff"]=new TProfile("eff", "eff", 400, 0, 20000);
  prof["eff_toppt"]=new TProfile("eff_toppt", "eff", 400, 0, 20000);
  prof["eff_leading_top"]=new TProfile("eff_leading_top", "eff", 400, 0, 20000);
  prof["eff_second_top"]=new TProfile("eff_second_top", "eff", 400, 0, 20000);

  prof["eff_leading_bottom"]=new TProfile("eff_leading_bottom", "eff", 400, 0, 20000);
  prof["eff_second_bottom"]=new TProfile("eff_second_bottom", "eff", 400, 0, 20000);

  prof["topfrac_toppt"]=new TProfile("topfrac_toppt", "top p_{T} fraction vs. top p_{T}", 400, 0, 20000);

  prof["toprecoilfrac_toppt"]=new TProfile("toprecoilfrac_toppt", "top recoil fraction vs. top p_{T}", 400, 0, 20000);

  prof["topfrac_jetpt"]=new TProfile("topfrac_jetpt", "top p_{T} fraction vs. jet p_{T}", 400, 0, 20000);

  prof["toprecoilfrac_jetpt"]=new TProfile("toprecoilfrac_jetpt", "top recoil fraction vs. jet p_{T}", 400, 0, 20000);

  prof["topfrac_met"]=new TProfile("topfrac_met", "top p_{T} fraction vs. MEt", 400, 0, 20000);

  prof["toppt_met"]=new TProfile("toppt_met", "top p_{T} vs. MEt", 400, 0, 20000);

  prof["eff_met"]=new TProfile("eff_met", "efficiency vs. MEt", 400, 0, 20000);

  prof["eff2_met"]=new TProfile("eff2_met", "efficiency vs. MEt (both jets)", 400, 0, 20000);


  prof["toprecoilfrac_met"]=new TProfile("toprecoilfrac_met", "top recoil fraction vs. MEt", 400, 0, 20000);
				     
  prof["topeta1_met"]=new TProfile("topeta1_met", "leading top #eta vs. Missing E_{T}", 400, 0, 20000);

  prof["topeta2_met"]=new TProfile("topeta2_met", "second top #eta vs. Missing E_{T}", 400, 0, 20000);

  prof["toppt1_met"]=new TProfile("toppt1_met", "leading top p_{T} vs. Missing E_{T}", 400, 0, 20000);

  prof["toppt2_met"]=new TProfile("toppt2_met", "second top p_{T} vs. Missing E_{T}", 400, 0, 20000);
  


  hist2D["topfrac_met"]=new TH2D("topfrac_v_met", "top p_{T} fraction vs. MEt", 400, 0, 20000, 400, 0, 1.0);

  hist2D["toprecoilfrac_met"]=new TH2D("toprecoilfrac_v_met", "top recoil fraction vs. MEt", 400, 0, 20000, 0, 2.0);

  //  hist["nj100_after"]=new TH1D("nj100_after", "number of jets > 100 GeV", 10, 0, 10);

  map<string, TH1D*>::iterator it=hist.begin();
  for(;it!=hist.end(); ++it)
    {
      it->second->Sumw2();
    }

  map<string, TH2D*>::iterator it2=hist2D.begin();
  for(;it2!=hist2D.end(); ++it2)
    {
      it2->second->Sumw2();
    }

  map<string, TProfile*>::iterator 
    ip=prof.begin();
  for(;ip!=prof.end(); ++ip)
    {
      ip->second->Sumw2();
    }




  //initialize tree reading
  MyTree c(input);

  double weight_scale=1.0;
  
  float min_pt=100.;

  //cut_flow numbers
  double cutflow[5]={0};

  
  cout<<"Processing "<<c.t->GetEntries()<<" events..."<<endl;
  int division=ceil(c.t->GetEntries()/20.);

  //define covariant matrix:
  Covariance co(4);
  
  //the fisher discriminant
  double fisher_w[4]={-0.160772, 0.0808241, 1.15165, -0.518156};

  /*
    {0.100031, 0.291078, 1.35676, -0.181478, 1.91652, -0.294906, -0.137145, -0.383801, -0.222116, 0.134492};
  */
  double fisher_offset= -0.0823686;

  for(Long64_t i=0; i< c.t->GetEntries(); i++)
    {
      if(i% division == 0)
	cout<<(i/division)*5<<" % completed"<<endl;

      if(debug)
	cout<<"at event: "<<i<<endl;

      c.t->GetEntry(i);

      double weight=c.weight;

      //if scale greater than zero
      //it means this is signal, and weight 
      //is computed from weight_scale
      if(c.weight_scale > 0)
	weight = c.weight_scale;

      //quick nmuons and Ws
      //and W and Z
      int nmu_gen=0;
      int ne_gen=0;
      int ntau_gen=0;
      int nneu_gen=0;
      int nW_gen=0;
      int nZ_gen=0;

      vector<TLorentzVector> vlep_gen;
      vector<TLorentzVector> vneu_gen;
      vector<TLorentzVector> vb_gen;
      vector<TLorentzVector> vW1lep_gen;
      vector<TLorentzVector> vW2lep_gen;

      for(int i=0; i<50; ++i)
	if(c.statusTruth[i] == 3)
	  {
	    TLorentzVector j;
	    j.SetPtEtaPhiM(c.pTTruth[i], c.etaTruth[i], 
			   c.phiTruth[i], c.massTruth[i]);

	    

	    if(abs(c.pdgIdTruth[i]) == 5)
	      {
		vb_gen.push_back(j);
	      }

	    if(abs(c.pdgIdTruth[i]) == 11 || 
	       abs(c.pdgIdTruth[i]) == 13 ||
	       abs(c.pdgIdTruth[i]) == 15 )
	      {
		vlep_gen.push_back(j);
	      }

	    if(abs(c.pdgIdTruth[i]) == 12 || 
	       abs(c.pdgIdTruth[i]) == 14 ||
	       abs(c.pdgIdTruth[i]) == 16 )
	      {
		vneu_gen.push_back(j);
	      }

	    if(abs(c.pdgIdTruth[i]) == 13)
	      nmu_gen ++;

	    if(abs(c.pdgIdTruth[i]) == 11)
	      ne_gen ++;

	    if(abs(c.pdgIdTruth[i]) == 15)
	      ntau_gen ++;

	    if(abs(c.pdgIdTruth[i]) == 12 ||
	       abs(c.pdgIdTruth[i]) == 14 ||
	       abs(c.pdgIdTruth[i]) == 16)
	      nneu_gen ++;
	  
	    if(abs(c.pdgIdTruth[i]) == 24)
	      nW_gen ++;

	    if(abs(c.pdgIdTruth[i]) == 23)
	      nZ_gen ++;

	  }

      //W/Z cut
      if(b_W_only && nW_gen ==0)
	continue;

      if(b_Z_only && nZ_gen ==0)
	continue;

      //get the leptons from W
      if(vneu_gen.size()==2 && vlep_gen.size()==2)
	{
	  double mass11=
	    fabs((vneu_gen[0]+vlep_gen[0]).M() - 80.385);
	  double mass21=
	    fabs((vneu_gen[1]+vlep_gen[1]).M() - 80.385);

	  double mass12=
	    fabs((vneu_gen[0]+vlep_gen[1]).M() - 80.385);
	  double mass22=
	    fabs((vneu_gen[1]+vlep_gen[0]).M() - 80.385);

	  double metric1= mass11*mass11 + mass21*mass21;
	  double metric2= mass21*mass21 + mass22*mass22;
	  
	  
	  if(metric1 < metric2)
	    {
	      vW1lep_gen.push_back(vneu_gen[0]);
	      vW1lep_gen.push_back(vlep_gen[0]);

	      vW2lep_gen.push_back(vneu_gen[1]);
	      vW2lep_gen.push_back(vlep_gen[1]);
	    }

	  else
	    {
	      vW1lep_gen.push_back(vneu_gen[0]);
	      vW1lep_gen.push_back(vlep_gen[1]);

	      vW2lep_gen.push_back(vneu_gen[1]);
	      vW2lep_gen.push_back(vlep_gen[0]);
	    }

	}



      hist["nmu_gen"]->Fill(nmu_gen, weight);
      hist["ne_gen"]->Fill(ne_gen, weight);
      hist["ntau_gen"]->Fill(ntau_gen, weight);
      hist["nneu_gen"]->Fill(nneu_gen, weight);
      hist["nW_gen"]->Fill(nW_gen, weight);
      hist["nZ_gen"]->Fill(nZ_gen, weight);

      

      //fast cut

      if(c.nJets100 < 2 || fabs(c.etaJet[0]) > 2.5 || 
	 fabs(c.etaJet[1])>2.5 )
	continue;


      /*
      if(c.nJets100 <2 || fabs(c.etaJet[0]) > 2.5 || 
	 fabs(c.etaJet[1])>2.5 || 
	 fabs(c.TrimmedMassJet[0]-175) > 50 ||
	 fabs(c.TrimmedMassJet[1]-160) > 65 ||
	 c.MET < 800. ||
	 fixphi(c.mindPhiMetJet100) < 1.0 ||
	 c.nLeps35 > 0
	 //c.METsignif < 40.
	 )
	continue;
      */

      //apply a muon cut
      vector<TLorentzVector> muons=
	c.get_muons(100.);

      if(muons.size()>0)
	hist["muon_pt1"]->Fill(muons[0].Pt(), weight);

      else
	hist["muon_pt1"]->Fill(-1.0, weight);
	
      
      vector<TLorentzVector> tops=
	c.get_particle(6);

      double sumpt_all=0;

      for(int ijet=0; ijet<c.nJets30; ijet++)
	sumpt_all+=c.pTJet[ijet];

      if(tops.size()>0)
	hist["top_pt1"]->Fill(tops[0].Pt(), weight);
      
      else
	hist["top_pt1"]->Fill(-1.0, weight);

	


      vector<TLorentzVector> bottoms=
	c.get_particle(5);

      //iterate over jets
      //and match them to muons
      int nmuj=0; 
      int nmuj_200=0;
      int nmuj_500=0;

      TLorentzVector mu1;
      double muj_ptratio=0;

      bool b_leadingjet_muon=false;
      bool b_leading_secondjet_muon=false;


      for(int m=0; m<muons.size(); ++m)
	for(int ijet=0; ijet<1; ijet++)
	if(c.pTJet[ijet]>500. && fabs(c.etaJet[ijet])<2.5)
	  {
	    TLorentzVector j;
	    j.SetPtEtaPhiM(c.pTJet[ijet], c.etaJet[ijet], 
			   c.phiJet[ijet], c.massJet[ijet]);

	    if(muons[m].DeltaR(j) > 0.5)
	      continue;
	    
	    if(muons[m].Pt() > 500.)
	      nmuj_500++;

	    if(muons[m].Pt() > 200)
	      nmuj_200++;

	    nmuj++;

	    if(ijet==0 && muons[m].Pt() > 200.)
	      b_leadingjet_muon=true;

	    if((ijet==0 || ijet==1) && muons[m].Pt() > 200.)
	      b_leading_secondjet_muon=true;
	    
	    if(nmuj==1)
	      {
		mu1=muons[m];
		muj_ptratio= mu1.Pt()/j.Pt();
	      }
	    break;
	  }
      
      prof["eff"]->Fill(c.pTJet[0] , b_leadingjet_muon,
			weight);


      if(tops.size()>0)
	prof["eff_toppt"]->Fill(tops[0].Pt() , b_leadingjet_muon,
				weight);


      if(tops.size()>1 && sumpt_all > 0 && 
	 c.pTJet[0] > 1000. && c.pTJet[1] > 1000.)
	{
	  hist["top_pt2"]->Fill(tops[1].Pt(), weight);
	  hist["top_dR"]->Fill(tops[0].DeltaR(tops[1]), weight);
	  prof["topfrac_toppt"]->Fill
	    (  tops[0].Pt() + tops[1].Pt(),
	       (tops[0].Pt() + tops[1].Pt())/sumpt_all,
	       weight );
	  prof["toprecoilfrac_toppt"]->Fill
	    (  tops[0].Pt() + tops[1].Pt(),
	       (tops[0] + tops[1]).Pt()/
	       (tops[0].Pt() + tops[1].Pt()),
	       weight );

	  prof["topfrac_met"]->Fill
	    (  c.MET,
	       (tops[0].Pt() + tops[1].Pt())/sumpt_all,
	       weight );
	  prof["toprecoilfrac_met"]->Fill
	    (  c.MET,
	       (tops[0] + tops[1]).Pt()/
	       (tops[0].Pt() + tops[1].Pt()),
	       weight );

	  prof["topfrac_jetpt"]->Fill
	    (  c.pTJet[0],
	       (tops[0].Pt() + tops[1].Pt())/sumpt_all,
	       weight );
	  prof["toprecoilfrac_jetpt"]->Fill
	    (  c.pTJet[0],
	       (tops[0] + tops[1]).Pt()/
	       (tops[0].Pt() + tops[1].Pt()),
	       weight );
	  
	  prof["toppt_met"]->Fill
	    (  c.MET,
	       fabs(tops[0].Pt()),
	       weight );

	  prof["eff_met"]->Fill
	    (  c.MET,
	       b_leadingjet_muon,
	       weight );

	  prof["eff2_met"]->Fill
	    (  c.MET,
	       b_leading_secondjet_muon,
	       weight );


	  prof["topeta1_met"]->Fill
	    (  c.MET,
	       fabs(tops[0].Eta()),
	       weight );

	  prof["topeta2_met"]->Fill
	    (  c.MET,
	       fabs(tops[1].Eta()),
	       weight );

	  prof["toppt1_met"]->Fill
	    (  c.MET,
	       tops[0].Pt(),
	       weight );

	  prof["toppt2_met"]->Fill
	    (  c.MET,
	       tops[1].Pt(),
	       weight );

	}

      else
	hist["top_pt2"]->Fill(-1.0, weight);



      //matching jets to top/bottom
      bool leading_bottom=false;
      bool second_bottom=false;
      bool leading_top=false;
      bool second_top=false;

	for(int ijet=0; ijet<2; ijet++)
	  if(c.pTJet[ijet]>100. && fabs(c.etaJet[ijet])<2.5)
	    {
	      TLorentzVector j;
	      j.SetPtEtaPhiM(c.pTJet[ijet], c.etaJet[ijet], 
			     c.phiJet[ijet], c.massJet[ijet]);
	      if(ijet==0)
		{
		  leading_bottom=match_jet(j, bottoms);
		  leading_top=match_jet(j, tops);
		}

	      if(ijet==1)
		{
		  second_bottom=match_jet(j, bottoms);
		  second_top=match_jet(j, tops);
		}

	    }

	prof["eff_leading_top"]->Fill
	  (c.pTJet[0], leading_top, weight);

	prof["eff_second_top"]->Fill
	  (c.pTJet[1], second_top, weight);


	prof["eff_leading_bottom"]->Fill
	  (c.pTJet[0], leading_bottom, weight);

	prof["eff_second_bottom"]->Fill
	  (c.pTJet[1], second_bottom, weight);
	


      //initialize variables
      int nj100=0;
      int nj150=0;
      //number of bjets
      int nbj100=0;

      int nl35=c.nLeps35;
      int nl100=0;
      int nl200=0;

      double HT=c.HTjets;
      double MJ=0;

      //get the met
      double met=c.MET;
      TLorentzVector vmet;
      vmet.SetPtEtaPhiM(met, 0, c.MET_phi, 0);
      
      double metsig=c.METsignif;
      double dphimin_500=999.;
      double dphimin_200=999.;
      double dphilep=999.;


      double dphi=fabs(fixphi(c.mindPhiMetJet100));

      //leading and second jet
      TLorentzVector j1, j2, j3, j4, l1, l2;
      double j1_t31, j1_t32, j1_tm, j1_dphi;
      double j2_t31, j2_t32, j2_tm, j2_dphi;

      //get input jets for MR computation
      vector<TLorentzVector> input_jets;
      
      int nl=0;

      //iterate over leptons
      for(int ilep=0; ilep<10; ilep++)
	if(c.pTLep[ilep]>=35.)
	  if(fabs(c.etaLep[ilep])<2.5)
	    {
	      TLorentzVector l;
	      l.SetPtEtaPhiM(c.pTLep[ilep], c.etaLep[ilep], 
			     c.phiLep[ilep], 0);
	      
	      
	      if(nl==0)
		l1=l;

	      else if(nl==1)
		l2=l;

	      if(fixphi(vmet.DeltaPhi(l)) < dphilep)
		dphilep=fixphi(vmet.DeltaPhi(l));

	      if(c.pTLep[ilep]>=100.)
		nl100++;

	      if(c.pTLep[ilep]>=200.)
		nl200++;
	      
	      nl++;
	    }      
      
      //iterate over jets
      for(int ijet=0; ijet<10; ijet++)
	if(c.pTJet[ijet]>=50.)
	  if(fabs(c.etaJet[ijet])<2.5)
	  {
	    if(debug)
	      cout<<"jet: "<<ijet<<endl;

	    TLorentzVector j;
	    j.SetPtEtaPhiM(c.pTJet[ijet], c.etaJet[ijet], 
			   c.phiJet[ijet], c.massJet[ijet]);
	    
	    input_jets.push_back(j);

	    //calcalate a dphi min
	    if(c.pTJet[ijet] >= 500)
	      if(fixphi(vmet.DeltaPhi(j)) < dphimin_500)
		dphimin_500=fixphi(vmet.DeltaPhi(j));

	    if(c.pTJet[ijet] >= 200)
	      if(fixphi(vmet.DeltaPhi(j)) < dphimin_200)
		dphimin_200=fixphi(vmet.DeltaPhi(j));
	      
	    if(c.pTJet[ijet] >= min_pt )//&& fabs(jet->Eta) < 3.0)
	      {
		nj100++;
		
		if(c.btagJet[ijet])
		  nbj100++;

		MJ+=c.massJet[ijet];
		
		//define substructure variables
		double c_tm=c.TrimmedMassJet[ijet];
		double c_t1=c.Tau1Jet[ijet];
		double c_t2=c.Tau2Jet[ijet];
		double c_t3=c.Tau3Jet[ijet];

		double c_t31=1;
		double c_t32=1;
		
		if(c_t1 > 0.01 && c_t2 > 0.01)
		  {
		    c_t31=c_t3/c_t1;
		    c_t32=c_t3/c_t2;
		  }

		/*
		double cphi=fabs(j.DeltaPhi(vmet));
		if(cphi < dphi)
		  dphi=cphi;
		*/

		switch(nj100)
		  {
		  case 1: 
		    j1=j;
		    j1_tm=c_tm;
		    j1_t32=c_t32;
		    j1_t31=c_t31;
		    break;
		      
		  case 2: 
		    j2=j;
		    j2_tm=c_tm;
		    j2_t32=c_t32;
		    j2_t31=c_t31;
		    break;
		    
		  case 3: 
		    hist["jet3pt"]->Fill(c.pTJet[ijet], weight);
		    j3=j;
		    break;
		    
		  case 4: 
		    hist["jet4pt"]->Fill(c.pTJet[ijet], weight);
		    j4=j;
		    break;
		  }
		
		if(c.pTJet[ijet] > 150)
		  nj150++;
	      }
	  }


      //compute dphi variables
      j1_dphi=fabs(vmet.DeltaPhi(j1));
      j2_dphi=fabs(vmet.DeltaPhi(j2));

      //change to the correct MR calculation
      vector<TLorentzVector> mega_jets=megajet_brute(input_jets);

      if(debug)
	cout<<"mega_jets size (should be 2): "
	    <<mega_jets.size()<<endl;

      double MR=cal_mr(mega_jets);

      //fill histograms
      hist["jet1pt"]->Fill(j1.Pt(), weight);
      hist["jet1mass"]->Fill(j1.M(), weight);
      hist["jet2pt"]->Fill(j2.Pt(), weight);
      hist["jet2mass"]->Fill(j2.M(), weight);


      hist["HT"]->Fill(HT, weight);
      hist["MJ"]->Fill(MJ, weight);
      hist["MR"]->Fill(MR, weight);
      hist["nj100"]->Fill(nj100, weight);
      hist["nbj100"]->Fill(nbj100, weight);
      hist["nj150"]->Fill(nj150, weight);
      hist["met"]->Fill(met, weight);
      hist["metsig"]->Fill(metsig, weight);
      hist["dphi"]->Fill(dphi, weight);
      hist["nmuj"]->Fill(nmuj, weight);


      //all but plots
      bool cut_pt=false;
      bool cut_mu=false;
      bool cut_mj1=true;
      bool cut_mj2=true;
      bool cut_met=false;
      bool cut_dphi=false;

      //only for high met_cut
      //the boosted regime
      if(met_cut > 10)
	{
	  if(c.pTJet[0] > 1000. && c.pTJet[1] > 1000.) cut_pt=true;
	  if(nmuj_200 > 0) cut_mu=true;
	  //if(fabs(c.TrimmedMassJet[0]-150) < 100.) cut_mj1=true;
	  //if(fabs(c.TrimmedMassJet[1]-150) < 100.) cut_mj2=true;
	  if(met > met_cut) cut_met=true;
	  if(dphimin_200 > 1.0 && dphimin_200 < 3.0) cut_dphi=true;
	}
      //now for compressed region
      else
	{
	  if(nl35 >= 2 && (l1+l2).M() < 120)
	    nl35=-1;

	  if(c.pTJet[0] > 200. && c.pTJet[1] > 200.) cut_pt=true;
	  if(nl35 >= 2) cut_mu=true;
	  //if(fabs(c.TrimmedMassJet[0]-150) < 100.) cut_mj1=true;
	  //if(fabs(c.TrimmedMassJet[1]-150) < 100.) cut_mj2=true;
	  if(met > 1500.) cut_met=true;
	  if(dphilep > .8 && dphimin_200 > 0.8) cut_dphi=true;
	}

      //allbut plots
      if(cut_mu && 
	 cut_met && cut_dphi)
	{
	  hist["jet1pt_allbut"]->Fill(c.pTJet[0], weight);
	  hist["jet2pt_allbut"]->Fill(c.pTJet[1], weight);
	  hist["HT_allbut"]->Fill(HT, weight);
	}

      if(cut_pt && 
	 cut_met && cut_dphi)
	{
	  hist["nmuj_allbut"]->Fill(nmuj, weight);
	  hist["nl35_allbut"]->Fill(nl35, weight);
	  hist["nl100_allbut"]->Fill(nl100, weight);
	  hist["nl200_allbut"]->Fill(nl200, weight);
	  hist["nmuj500_allbut"]->Fill(nmuj_500, weight);
	  hist["nmuj200_allbut"]->Fill(nmuj_200, weight);
	}

      if(cut_pt && cut_mu &&
	 cut_met && cut_dphi)
	hist["jet1masstrim_allbut"]->
	  Fill(c.TrimmedMassJet[0], weight);

      if(cut_pt && cut_mu &&
	 cut_met && cut_dphi)
	hist["jet2masstrim_allbut"]->
	  Fill(c.TrimmedMassJet[1], weight);


      if(cut_pt && 
	 cut_dphi && tops.size()>1)
	{
	  hist2D["topfrac_met"]->Fill	    
	    (  c.MET,
	       (tops[0].Pt() + tops[1].Pt())/sumpt_all,
	       weight );

	  hist2D["toprecoilfrac_met"]->Fill
	    (  c.MET,
	       (tops[0] + tops[1]).Pt()/
	       (tops[0].Pt() + tops[1].Pt()),
	       weight );
	} 
   
      if(cut_pt && cut_mu && 
	 cut_dphi && (c.nLeps35 == 0 || met_cut <=10) )
	{
	  hist["met_allbut"]->Fill(met, weight);
	  hist["metsig_allbut"]->Fill(metsig, weight);
	}


      if(cut_pt && nneu_gen==2 && vlep_gen.size()==2 &&
	 vb_gen.size()==2)
	{
	  //loop over all possible changes neu -> lep...etc
	  TLorentzVector met_mod=vmet-vneu_gen[0]-vneu_gen[1]
	    + vlep_gen[0] + vlep_gen[1];

	  for(int l1=0; l1<2; ++l1)
	    for(int l2=0; l2<2; ++l2)
	      {
		int n1=(l1+1) %2;
		int n2=(l2+1) %2;

		TLorentzVector met_new=
		  met_mod - vW1lep_gen[l1] - vW2lep_gen[l2]
		  + vW1lep_gen[n1] + vW2lep_gen[n2];
		
		//get the new deltaphi
		double dphimin_mod=999;
		for(int ijet=0; ijet<10; ijet++)
		  if(c.pTJet[ijet] >= 200 && 
		     fabs(c.etaJet[ijet])<2.5)
		    {
		      TLorentzVector j;
		      j.SetPtEtaPhiM(c.pTJet[ijet], 
				     c.etaJet[ijet], 
				     c.phiJet[ijet], 
				     c.massJet[ijet]);

		      if(fixphi(met_new.DeltaPhi(j)) 
			 < dphimin_mod)
			dphimin_mod=fixphi(met_new.DeltaPhi(j));
		    }

		//get the new deltaphilep
		double dphilep_mod1=
		  fixphi(met_new.DeltaPhi(vW1lep_gen[l1]));

		double dphilep_mod2=
		  fixphi(met_new.DeltaPhi(vW2lep_gen[l2]));
		
		if(dphilep_mod2 < dphilep_mod1)
		  dphilep_mod1 = dphilep_mod2;

		double min_dR=999;

		  for(int k=0; k<2; ++k)
		    {
		      double temp_dR1=
			vW1lep_gen[l1].DeltaR(vb_gen[k]);
		      double temp_dR2=
			vW2lep_gen[l2].DeltaR(vb_gen[k]);
		      
		      if(temp_dR1 < min_dR)
			min_dR = temp_dR1;

		      if(temp_dR2 < min_dR)
			min_dR = temp_dR2;
		    }



		if(min_dR > 0.5 && min_dR < 999 &&
		   (vW1lep_gen[l1]+vW2lep_gen[l2]).M() > 120 &&
		   dphilep_mod1 > 0.8 &&
		   dphimin_mod > 0.8 &&		   
		   vW1lep_gen[l1].Pt() > 35 &&
		   vW2lep_gen[l2].Pt() > 35 &&
		   fabs(vW1lep_gen[l1].Eta()) < 2.5 &&
		   fabs(vW2lep_gen[l2].Eta()) < 2.5
		   )
		  hist["met_allbut_dilep_mod"]->Fill(met_new.Pt(), weight);

		
		
	      }
	    
	  

	}


      if(nneu_gen==2 && cut_pt)
	{
	  if(vlep_gen.size()==2 && 
	     (vlep_gen[0]+vlep_gen[1]).M() > 120 &&
	     fixphi(vmet.DeltaPhi(vlep_gen[0])) > 0.8 &&
	     fixphi(vmet.DeltaPhi(vlep_gen[1])) > 0.8 && 
	     vlep_gen[0].Pt() > 35 &&
	     vlep_gen[1].Pt() > 35 &&
	     fabs(vlep_gen[0].Eta()) < 2.5 &&
	     fabs(vlep_gen[1].Eta()) < 2.5)
	    hist["met_allbut_dilep_nolep"]->Fill(met, weight);
	  
	  double min_dR=999;
	  if(vlep_gen.size()==2 && vb_gen.size()==2)
	    for(int k=0; k<2; ++k)
	      for(int l=0; l<2; ++l)
		{
		  double temp_dR=vlep_gen[k].DeltaR(vb_gen[l]);
		  if(temp_dR < min_dR)
		    min_dR = temp_dR;
		}

	  if(min_dR > 0.5 && min_dR < 999 &&
	     (vlep_gen[0]+vlep_gen[1]).M() > 120 &&
	     fixphi(vmet.DeltaPhi(vlep_gen[0])) > 0.8 && 
	     fixphi(vmet.DeltaPhi(vlep_gen[1])) > 0.8 && 
	     vlep_gen[0].Pt() > 35 &&
	     vlep_gen[1].Pt() > 35 &&
	     fabs(vlep_gen[0].Eta()) < 2.5 &&
	     fabs(vlep_gen[1].Eta()) < 2.5
	     )
	    {
	      if(dphimin_200 > 0.8)
		hist["met_allbut_dilep_dR"]->Fill(met, weight);


	      if(dphimin_200 > 0.5)
		hist["met_allbut_dilep_dR_loose"]->Fill(met, weight);
	    }

	}

      if(cut_pt && cut_mu && 
	 cut_dphi)
	{
	  if(nneu_gen==2)
	    hist["met_allbut_dilep"]->Fill(met, weight);


	  hist["nmu_after"]->Fill(nmu_gen, weight);
	  hist["ne_after"]->Fill(ne_gen, weight);
	  hist["ntau_after"]->Fill(ntau_gen, weight);
	  hist["nneu_after"]->Fill(nneu_gen, weight);
	}

      if(cut_pt && 
	 cut_dphi )
	{
	  hist["met_allbut_nolep"]->Fill(met, weight);
	}

      if(cut_pt && cut_mu && dphimin_200 > 0.5 )
	{
	  hist["met_allbut_mild"]->Fill(met, weight);
	}

      if(cut_pt && dphimin_200 > 0.3)
	{
	  hist["met_allbut_milder"]->Fill(met, weight);
	}

      if(cut_pt && cut_mu && cut_mj1 && cut_mj2 &&
	 cut_met)
	{
	  hist["dphimin200_allbut"]->Fill(dphimin_200, weight);
	  hist["dphimin500_allbut"]->Fill(dphimin_500, weight); 
	  hist["dphilep_allbut"]->Fill(dphilep, weight); 
	}

      //basic cutflow
      hist["cutflow"]->Fill(0., weight);
      if(cut_pt){
	hist["cutflow"]->Fill(1., weight);

      if(cut_mu){
	hist["cutflow"]->Fill(2., weight);
	
	/*
	//175, 50
	if(cut_mj1){
	hist["cutflow"]->Fill(3., weight);
	
	//160, 65
	
	if(cut_mj2){
	  hist["cutflow"]->Fill(4., weight);
	*/
  
	//if(cut_met){
	if(c.nLeps35 == 0 || met_cut < 10){
	    hist["cutflow"]->Fill(3., weight);
	
	    //if(cut_dphi){
	    if(cut_dphi){
	      hist["cutflow"]->Fill(4., weight);

	      //if(c.nLeps35 == 0 || met_cut < 10){
	      if(cut_met){
		hist["cutflow"]->Fill(5., weight);


		hist["HT_after"]->Fill(HT, weight);
		hist["MR_after"]->Fill(MR, weight);
		hist["nj100_after"]->Fill(nj100, weight);
		hist["nbj100_after"]->Fill(nbj100, weight);
		hist["met_after"]->Fill(met, weight);
		hist["metsig_after"]->Fill(metsig, weight);
		hist["dphi_after"]->Fill(dphi, weight);
		hist["jet1pt_after"]->Fill(j1.Pt(), weight);
		hist["jet2pt_after"]->Fill(j2.Pt(), weight);
		hist["jet1mass_after"]->Fill(j1.M(), weight);
		hist["jet2mass_after"]->Fill(j2.M(), weight);	
		hist["mjj_after"]->Fill((j1+j2).M(), weight);
		hist["mll_after"]->Fill((l1+l2).M(), weight);
		
		hist["ptl1_after"]->Fill(l1.Pt(), weight);
		hist["ptl2_after"]->Fill(l2.Pt(), weight);	

		hist["dphill_after"]->Fill
		  (fixphi(l1.DeltaPhi(l2)), weight);	

		hist["dRll_after"]->Fill
		  (l1.DeltaR(l2), weight);	

		hist["ST_after"]->Fill
		  (l1.Pt() + l2.Pt() + 
		   j1.Pt() + j2.Pt() + j3.Pt() +
		   met, weight);

		hist["nj100_after"]->Fill
		  (nj100, weight);

		if(tops.size()>0)
		  hist["top_pt1_after"]->Fill(tops[0].Pt(), weight);

		if(tops.size()>1)
		  {
		    hist["top_dR_after"]->Fill
		      (tops[0].DeltaR(tops[1]), weight);
		    hist["top_pt2_after"]->Fill(tops[1].Pt(), weight);
		  }

		
	      }}}}}


      //fill cutflow
      hist["cutflow2"]->Fill(0., weight);
      if(c.pTJet[0] > 1000.){
	hist["cutflow2"]->Fill(1., weight);

      if(nmuj > 0){
	hist["cutflow2"]->Fill(2., weight);
	
	//175, 50
	//if(fabs(c.TrimmedMassJet[0]-150) < 100.){
	if(true){
	hist["cutflow2"]->Fill(3., weight);
	
	//160, 65
	//if(fabs(c.TrimmedMassJet[1]-150) < 100.){
	if(true){
	  hist["cutflow2"]->Fill(4., weight);
	  
	  if(met > 3000.){
	    hist["cutflow2"]->Fill(5., weight);

	    if(dphimin_200 > 1.0){
	      hist["cutflow2"]->Fill(6., weight);

	      if(c.nLeps35 == 0){
		hist["cutflow2"]->Fill(7., weight);

	      }}}}}}}
      

      //fill cutflow
      hist["cutflow3"]->Fill(0., weight);
      if(c.pTJet[0] > 1000.){
	hist["cutflow3"]->Fill(1., weight);

      if(nmuj > 0){
	hist["cutflow3"]->Fill(2., weight);
	
	//175, 50
	//if(fabs(c.TrimmedMassJet[0]-150) < 100.){
	if(true){
	hist["cutflow3"]->Fill(3., weight);
	
	//160, 65
	//if(fabs(c.TrimmedMassJet[1]-150) < 100.){
	if(true){
	  hist["cutflow3"]->Fill(4., weight);
	  
	  if(metsig > 45.){
	    hist["cutflow3"]->Fill(5., weight);

	    if(dphimin_500 > 1.0){
	      hist["cutflow3"]->Fill(6., weight);

	      if(c.nLeps35 == 0){
		hist["cutflow3"]->Fill(7., weight);

	      }}}}}}}

      //fill cutflow
      hist["cutflow4"]->Fill(0., weight);
      if(c.pTJet[0] > 500.){
	hist["cutflow4"]->Fill(1., weight);

      if(nmuj >= 2){
	hist["cutflow4"]->Fill(2., weight);
	
	//175, 50
	//if(fabs(c.TrimmedMassJet[0]-150) < 100.){
	if(true){
	hist["cutflow4"]->Fill(3., weight);
	
	//160, 65
	//if(fabs(c.TrimmedMassJet[1]-150) < 100.){
	if(true){
	  hist["cutflow4"]->Fill(4., weight);
	  
	  if(metsig > 45.){
	    hist["cutflow4"]->Fill(5., weight);

	    if(dphimin_200 > 1.0){
	      hist["cutflow4"]->Fill(6., weight);

	      if(c.nLeps35 == 0){
		hist["cutflow4"]->Fill(7., weight);

	      }}}}}}}
      

      //fill cutflow
      hist["cutflow5"]->Fill(0., weight);
      if(c.pTJet[0] > 500.){
	hist["cutflow5"]->Fill(1., weight);

      if(nmuj > 1){
	hist["cutflow5"]->Fill(2., weight);
	
	//175, 50
	//if(fabs(c.TrimmedMassJet[0]-150) < 100.){
	if(true){
	hist["cutflow5"]->Fill(3., weight);
	
	//160, 65
	//if(fabs(c.TrimmedMassJet[1]-150) < 100.){
	if(true){
	  hist["cutflow5"]->Fill(4., weight);
	  
	  if(metsig > 45.){
	    hist["cutflow5"]->Fill(5., weight);

	    if(dphimin_200 > 1.0){
	      hist["cutflow5"]->Fill(6., weight);

	      if(c.nLeps35 == 0){
		hist["cutflow5"]->Fill(7., weight);

	      }}}}}}}


      hist["cutflow6"]->Fill(0., weight);
      if(c.pTJet[0] > 1000.){
	hist["cutflow6"]->Fill(1., weight);

      if(nmuj_500 > 0){
	hist["cutflow6"]->Fill(2., weight);
	
	//175, 50
	//if(fabs(c.TrimmedMassJet[0]-150) < 100.){
	if(true){
	hist["cutflow6"]->Fill(3., weight);
	
	//160, 65
	//if(fabs(c.TrimmedMassJet[1]-150) < 100.){
	if(true){
	  hist["cutflow6"]->Fill(4., weight);
	  
	  if(metsig > 45.){
	    hist["cutflow6"]->Fill(5., weight);

	    if(dphimin_200 > 1.0){
	      hist["cutflow6"]->Fill(6., weight);

	      if(c.nLeps35 == 0){
		hist["cutflow6"]->Fill(7., weight);

	      }}}}}}}
      

      



      //after some mild cut
      //MR > 1000
      //met > 500
      //if(MR > 1000 && met > 500)
      if(c.nJets100 >=2 && 
	 fabs(c.etaJet[0]) < 2.5 &&
	 fabs(c.etaJet[1]) < 2.5 &&
	 //fabs(c.TrimmedMassJet[0]-175) < 75 &&
	 //fabs(c.TrimmedMassJet[1]-175) < 75 &&
	 //c.MET > 600. &&
	 dphimin_200 > 1.0
	 //c.nLeps35 == 0 &&
	 //nmuj_500 > 0
	 )
	{
	  hist["HT_mild"]->Fill(HT, weight);
	  hist["MR_mild"]->Fill(MR, weight);
	  hist["MRHT_mild"]->Fill(MR/HT, weight);
	  hist["nbj100_mild"]->Fill(nbj100, weight);
	  
	  hist["met_mild"]->Fill(met, weight);
	  hist["metsig_mild"]->Fill(metsig, weight);
	  
	  hist["dphi1_mild"]->Fill(fixphi(j1_dphi), weight);
	  hist["dphi2_mild"]->Fill(fixphi(j2_dphi), weight);

	  hist["dphimin_mild"]->Fill(fixphi(c.mindPhiMetJet100), weight);

	  hist["dphimin500_mild"]->Fill(fixphi(dphimin_500), weight);
	  hist["dphimin500_mod_mild"]->Fill(transformphi(dphimin_500), weight);

	  hist["dphimin200_mild"]->Fill(fixphi(dphimin_200), weight);
	  hist["dphimin200_mod_mild"]->Fill(transformphi(dphimin_200), weight);

	  hist["dphi1_mod_mild"]->Fill(transformphi(j1_dphi), weight);
	  hist["dphi2_mod_mild"]->Fill(transformphi(j2_dphi), weight);

	  hist["dphimin_mod_mild"]->Fill
	    (transformphi(c.mindPhiMetJet100), weight);
	  
	  hist["dR_mild"]->Fill(j1.DeltaR(j2), weight);
	  hist["dphij_mild"]->Fill(fabs(j1.DeltaPhi(j2)), weight);
	  hist["deta_mild"]->Fill(fabs(j1.Eta()-j2.Eta()), weight);

	  hist["ptmuj_mild"]->Fill(muj_ptratio, weight);

	  hist["jetptratio_mild"]->Fill(j2.Pt()/j1.Pt(), weight);

	  hist["jet1pt_mild"]->Fill(j1.Pt(), weight);
	  hist["jet1mass_mild"]->Fill(j1.M(), weight);
	  hist["jet1masstrim_mild"]->Fill(j1_tm, weight);
	  hist["jet1t32_mild"]->Fill(j1_t32, weight);
	  hist["jet1t31_mild"]->Fill(j1_t31, weight);

	  hist["jet2pt_mild"]->Fill(j2.Pt(), weight);
	  hist["jet2mass_mild"]->Fill(j2.M(), weight);
	  hist["jet2masstrim_mild"]->Fill(j2_tm, weight);
	  hist["jet2t32_mild"]->Fill(j2_t32, weight);
	  hist["jet2t31_mild"]->Fill(j2_t31, weight);

	  hist["nj100_mild"]->Fill(nj100, weight);
	  hist["nl35_mild"]->Fill(c.nLeps35, weight);
	  hist["mjj_mild"]->Fill((j1+j2).M(), weight);
	  
	  double pt_norm=c.pTJet[0];
	  
	  double var_list[]=
	    {
	      log(c.pTJet[0]/min_pt), 
	      log(HT/500.),
	      //log(1-c.pTJet[1]/c.pTJet[0]),
	      log(met/1000.),
	      log(metsig/40)
	      //j1_t32,
	      //j2_t32
	    };
	  
	  int list_size= sizeof(var_list)/sizeof(double);
	  assert(list_size==co.Rank());
	  
	  vector<double> list(var_list, 
			      var_list+ list_size);
	  
	  co.Fill(list, weight);
	  
	  //compute the fisher discriminant
	  double fisher=0;
	  for(int k=0; k<list_size; ++k)
	    {
	      fisher+=list[k]*fisher_w[k];
	    }
	  fisher+= fisher_offset;
	  
	  hist["fisher"]->Fill(fisher, weight);
	}
    

      //fill old cutflow plots      
      hist["cutflow_old"]->Fill(0., weight);
      if( c.massJet[0] <= 300 && c.massJet[0] >=100 ){
	hist["cutflow_old"]->Fill(1, weight);
	
      if( c.massJet[1] <= 300 && c.massJet[1] >=100 ){
	hist["cutflow_old"]->Fill(2, weight);
	  
      if(fabs(vmet.DeltaPhi(j1)) >= 1.6 && 
	 fabs(vmet.DeltaPhi(j2)) >= 0.8 ){
	    
	hist["cutflow_old"]->Fill(3, weight);
	    
      if(MR > MR_cut){
	hist["cutflow_old"]->Fill(4, weight);
		  
      if(met >= met_cut){
	hist["cutflow_old"]->Fill(5, weight);
	
      if(metsig >= 50){
	  hist["cutflow_old"]->Fill(6, weight);
	  
      if(nmuj >= 1){
	hist["cutflow_old"]->Fill(7, weight);		
      
      if(nmuj >= 2){
	hist["cutflow_old"]->Fill(8, weight);		
    }}}}}}}}
	

      
      /*
      //cuts
      if(nj100 >= 2 &&
	 nbj100 >= 2 &&
	 j2.M() >= 160 &&
	 dphi >= 0.32 &&
	 met >= 5700 &&
	 metsig >= 66)
	{
	  hist["HT_after"]->Fill(HT, weight);
	  
	  hist["nj100_after"]->Fill(nj100, weight);
	  hist["nbj100_after"]->Fill(nbj100, weight);
	  
	  hist["met_after"]->Fill(met, weight);
	  hist["metsig_after"]->Fill(metsig, weight);
	  hist["dphi_after"]->Fill(dphi, weight);

	  hist["jet1pt_after"]->Fill(j1.Pt(), weight);	  
	  hist["jet2pt_after"]->Fill(j2.Pt(), weight);	  

	  hist["jet1mass_after"]->Fill(j1.M(), weight);	  
	  hist["jet2mass_after"]->Fill(j2.M(), weight);	  
	}
      */

      //here do a fisher discriminant
      //the list of variables are:

      //pt1, pt2, HT, MR, met, met_sig, dphi1, dphi2,
      //t32_firstjet, t31_firstjet, dR between jets

    }

  cout<<"100 % completed"<<endl;

  cout<<co<<endl;


  cout<<"cut flow: ";

  int cutflowsize= sizeof(cutflow)/sizeof(double);
  for(int i=0; i< cutflowsize; ++i)
    {
      cout<<cutflow[i];

      if(i== cutflowsize-1)
	{
	  cout<<endl;
	  break;
	}
      
      cout<<", ";
    }

  TFile fout(output.c_str(), "RECREATE");
  //  map<string, TH1D*>::iterator 
  it=hist.begin();
  for(;it!=hist.end(); ++it)
    {
      it->second->Write();
    }

  it2=hist2D.begin();
  for(;it2!=hist2D.end(); ++it2)
    {
      it2->second->Write();
    }
  
  ip=prof.begin();
  for(;ip!=prof.end(); ++ip)
    {
      ip->second->Write();
    }

  fout.Close();
  return 1;
}
